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We investigate the emergence of singular solutions in a non-local model for a 
magnetic system. We study a modified Gilbert-type equation for the magnetization 
vector and find that the evolution depends strongly on the length scales of the non- 
local effects. We pass to a coupled density-magnetization model and perform a linear 
stability analysis, noting the effect of the length scales of non-locality on the system's 
stability properties. We carry out numerical simulations of the coupled system and 
find that singular solutions emerge from smooth initial data. The singular solutions 
represent a collection of interacting particles (clumpons). By restricting ourselves to 
the two-clumpon case, we are reduced to a two-dimensional dynamical system that 
is readily analyzed, and thus we classify the different clumpon interactions possible. 

I. INTRODUCTION 

In recent years, the modeling of nanoscale physics has become important, both because 
of industrial applications [H El EJ H] , and because of the development of experiments that 
probe these small scales [5]. One particular problem is the modeling of aggregation, in 
which microscopic particles collapse under the potential they exert on each other, and form 
mesoscopic structures that in turn behave like particles. 

In a series of papers, Holm, Putkaradze and Tronci |Hl El El El Ell ED] have focused on 
the derivation of aggregation equations that possess emergent singular solutions. Contin- 
uum aggregation equations have been used to model gravitational collapse and the subse- 
quent emergence of stars [12J, the localization of biological populations [131 HH US]; an d 
the self-assembly of nanoparticles [IB]. These are complexes of atoms or molecules that 
form mesoscale structures with particle-like behavior. The utility of the Holm-Putkaradze 
model lies in its emphasis on non-local physics, and the emergence of singular solutions from 
smooth initial data. Because of the singular (delta-function) behavior of the model, it is 
an appropriate way to describe the universal phenomena of aggregation and the subsequent 
formation of particle-like structures. Indeed in this framework, it is possible to prescribe the 
dynamics of the particle-like structures after collapse. Thus, the model provides a descrip- 
tion of directed self-assembly in nanophysics [TBI E], in which the detailed physics is less 
important than the effective medium properties of the dynamics. 

In this work we focus on equations introduced by Holm, Putkaradze and Tronci for 
the aggregation of oriented particles [Bl E]. We treat the initial state of the system as a 
continuum, a good approximation in nanophysics applications [18] . One realization of this 
problem is in nanomagnetics, in which particles with a definite magnetic moment collapse 
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and form mesoscale structures, that in turn have a definite magnetic moment. Thus, in this 
paper we refer to the orientation vector in our continuum picture as the magnetization. We 
investigate these equations numerically and study their evolution and aggregation properties. 
One aspect of non-local problems, already mentioned in [TTJ, is the effect of competition 
between the length scales of non-locality on the system evolution. We shall highlight this 
effect with a linear stability analysis of the full density-magnetization equations. 

This paper is organized as follows. In Sec. |n] we introduce a non-local Gilbert (NG) 
equation to describe non-local interactions in a magnetic system. We investigate the com- 
petition between the system's two length scales of non-locality. In Sec. Ill we introduce a 
coupled density-magnetization system that generates singular solutions. We examine the 
competition of length scales through a linear stability analysis and through the study of 
the dynamical equations for a simple singular solution that describes the interaction of two 
particle-like objects (clumpons). We perform numerical simulations that highlight the emer- 



gence of singular solutions from smooth initial data. We draw our conclusions in Sec. IV 



II. THE NON-LOCAL GILBERT EQUATION 

In this section we study a magnetization equation that in form is similar to the Gilbert 
equation, that is, the Landau-Lifshitz-Gilbert equation in the over-damped limit [THl EH] ■ 
The equation we focus on incorporates non-local effects, and was introduced in [B] . We study 
the evolution and energetics of this equation, and examine the importance of the problem 
length scales in determining the evolution. 

We study the following non-local Gilbert (NG) equation, 

d ( 5E\ 

m x x *zr > (!) 



dt \ 8m J 

where m is the magnetization density, fi m is the mobility, defined as 

^=(l-^) _1 m, 
and SE/Sm is the variational derivative of the energy, 

^=(l-o 2 ^r 1 m. 
5m v ' 

The smoothened magnetization \i m and the force d~E/5m can be computed using the theory 
of Green's functions. In particular, 



fi m (x,t)= / dyHp(x-y)m(y,t) := Hp * m(x,t) . 
Jn 

Here * denotes the convolution of functions, and the kernel Hp (x) satisfies the equation 

-f3 2 ^jHp(x) = 6(x). (2) 

The function 5 (x) is the Dirac delta function. Equation (J2| is solved subject to conditions 
imposed on the boundary of the domain Q. In this paper we shall work with a periodic 
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domain Q = [—L/2, L/2] or Q = [0, L], although other boundary conditions are possible. 
Note that Eq. has a family of non-trivial equilibrium states given by 



m eq (£J 



mo sin (fee + O ) 



where mo is a constant vector, k is some wave number, and 0o is a constant phase. The 
derivation of this solution is subject to the boundary conditions discussed in Sec. |III[ 

By setting j3 = and replacing (1 — a 2 d 2 ) with —d 2 , we recover the more familiar 
Landau-Lifshitz-Gilbert equation, in the overdamped limit 



dm 
~dt 



-m x m x 



d 2 m 
dx 2 



(3) 



Equation Q possesses several features that will be useful in understanding the numerical 
simulations. There is an energy functional 



E{t) = \ [ dxm ■ (1 - a 2 d 2 x ) ' 
Jn 



m, 



(4) 



which evolves in time according to the relation 



dE 
~dT 
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2n2\-! 



/i m • (1 - a 2 d 2 x ) 
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(l - a 2 dl) 1 m ■ fi m x(l-a 2 dl) 1 



(5) 



This is not necessarily a non-increasing function of time, although setting (3 = gives 



dE\ 
~dt) 



dx 

dxm 2 



m ■ (l — a 2 dl) 



dxm z 



n 2 



7» 



m 



(cos 2 — 1) < 0, 



(6) 

where ip is the angle between m and (1 — a 2 d 2 ) 1 m. In the special case when /9 — > 0, 
we therefore expect 75 (t) to be a non-increasing function of time. On the other hand, 
inspection of Eq. ^ shows that as a — ► 0, the energy tends to a constant. Additionally, the 
magnitude of the vector m is conserved. This can be shown by multiplying Eq. (fij) by m, 
and by exploiting the antisymmetry of the cross product. Thus, we are interested only in 
the orientation of the vector m; this can be parametrized by two angles on the sphere: the 
azimuthal angle 9 (x, t), and the polar angle (x, t), where 



I m I cosc^sin^, 



1 777 1 sin sin 6, 



m 7 _ 



\m\ cos6, 



(7) 



and where G [0, 2ir), and 9 G [0, n]. 

We carry out numerical simulations of Eqs. ([[J and ^ on a periodic domain [0,L], and 
outline the findings in what follows. Motivated by the change of coordinates (J7|, we choose 
the initial data 



0o (x) = 7r (1 + sin (2r7Tx/L)) , ^ (x) = |7r (1 + sin {2irsx/L)) , 



(8) 
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FIG. 1: (Color online) The initial data for the magnetization equation Q. This initialization 
is obtained by allowing the orientation angles of the magnetization vector to vary sinusoidally in 
space, as in Eq. Q. Here the wave number of the variation is equal to the fundamental wave 
number 2ir/L. 



where r and s are integers. These data are shown in Fig. [T[ 

Case 1: Numerical simulations of Eq. Equation Q is usually solved by explicit or 
implicit finite differences [20] . We solve the equation by these methods, and by the explicit 
spectral method [2T]. The accuracy and computational cost is roughly the same in each case, 
and for simplicity, we therefore employ explicit finite differences; it is this method we use 
throughout the paper. Given the initial conditions QSJ) , each component of the magnetization 
m = (m x ,m y ,m z ) tends to a constant, the energy 



E 



\ dx 


dm 


In 


dx 



decays with time, and \m\ retains its initial value \m\ 2 = 1. After some transience, the 
decay of the energy functional becomes exponential in time. These results are shown in 
Fig.|| 
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FIG. 2: (Color online) Numerical simulations of Case (1), the Landau-Lifshitz-Gilbert equation 
in the over-damped limit. In this case, the magnetization decays to a constant state. Subfigures 
(a) and (b) show the magnetization at times t = 0.03 and i = 0.15 respectively; (c) is the energy 
functional, which exhibits exponential decay after some transience. The final orientation is (</>, 9) = 
(7r,7r/2). 



Case 2: Numerical simulations of Eq. Q with a < (3. Given the smooth initial data (j8l), 
in time each component of the magnetization m = (m x ,m y ,m z ) decays to zero, while the 
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energy 

Jn 

tends to a constant value. Given our choice of initial conditions, the energy in fact increases 
to attain this constant value. Again the quantity |m| 2 stays constant. These results are 
shown in Fig. [3j We find similar results when we set a = 0. 



(a) 



0.2 0.4 0.6 0.8 



0.5 



_ _ _ m 

m 



-0.5 



(b) 




10 20 30 40 



(c) 



FIG. 3: (Color online) Numerical simulations of Case (2), the non-local Gilbert equation with with 
a < (3. In this case, the energy increases to a constant value, and the magnetization becomes 
constant. Subfigures (a) and (b) show the magnetization at times t = 8 and t = 40; (c) is the 
energy functional. The final orientation is ((f), 0) = (ir, 7r/2). 



Case 3: Numerical simulations of Eq. (JTJ with a > f3. Given the smooth initial data (J8J), 
in time each component of the magnetization m = (m x ,m y ,m z ) develops finer and finer 
scales. The development of small scales is driven by the decreasing nature of the energy 
functional, which decreases as power law at late times, and is reflected in snapshots of the 
power spectrum of the magnetization vector, shown in Fig. |4j As the system evolves, there 
is a transfer of large amplitudes to higher wave numbers. This transfer slows down at late 



Case 


Length scales 


Energy 


Outcome as t — > oo 


Linear Stability 


(1) 


P 


= 0, SE/5m = —d^m 


Decreasing 


Constant state 


Stable 


(2) 




a < (3 


Increasing 


Constant state 


Stable 


(3) 




a> [3 


Decreasing 


Development of finer and finer scales 


Unstable 



TABLE I: Summary of the forms of Eq. ([TJ studied. 



times, suggesting that the rate at which the solution roughens tends to zero, as t — > oo. 
The evolution preserves the symmetry of the magnetization vector m (x, t) under parity 
transformations. This is seen by comparing Figs. [T] and |4} The energy is a decaying function 
of time, while the quantity \m\ 2 stays constant. We find similar results for the case when 
(3 = 0. 

These results can be explained qualitatively as follows. In Case (1), the energy functional 
exacts a penalty for the formation of gradients. The energy decreases with time and the 
the system evolves into a state in which no magnetization gradients are present, that is, a 
constant state. On the other hand, we have demonstrated that in Case (2), when a < (3, 
the energy increases to a constant value. Since in the nondocal model, the energy functional 
represents the cost of forming smooth spatial structures, an increase in energy produces a 
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FIG. 4: (Color online) Numerical simulations of Case (3), the non-local Gilbert equation with 
a > (3. In this case, the energy decreases indefinitely, and the magnetization vector develops finer 
and finer scales. Subfigures (a), (b), and (c) show the magnetization at time t = 10000; (d) is the 
energy functional, which decreases in time as a power law at late times. Subfigure (e) shows the 
power spectrum of m x ; the integer index n labels the spatial scales: if k n is a wavenumber, then 
the corresponding integer label is n = k n L/2ir. 



smoother magnetization field, a process that continues until the magnetization reaches a 
constant value. Finally, in Case (3), when a > /3, the energy functional decreases, and this 
decrease corresponds to a roughening of the magnetization field, as seen in Fig. |4j In Sec. Ill 



we shall show that Case (2) is stable to small perturbations around a constant state, while 
Case (3) is unstable. Furthermore, we note that Case (2) and Case (3) differ only by a minus 
sign in Eq. Q, and are therefore related by time reversal. These results are summarized in 
Table E 

The solutions of Eqs. ([I]) and (|3) do not become singular. This is not surprising: the man- 
ifest conservation of |m| 2 in Eqs. (TJ and ^ provides a pointwise bound on the magnitude 
of the solution, preventing blow-up. Any addition to Eq. ([I]) that breaks this conservation 
law gives rise to the possibility of singular solutions, and it is to this possility that we now 
turn. 



III. COUPLED DENSITY-MAGNETIZATION EQUATIONS 

In this section we study a coupled density-magnetization equation pair that admit singular 
solutions. We investigate the linear stability of the equations and examine the conditions 
for instability. We find that the stability or otherwise of a constant state is controlled by the 
magnetization and density values of that state, and by the relative magnitude of the problem 
length scales. Using numerical and analytical techniques, we investigate the emergence and 
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self-interaction of singular solutions. 
The equations we study follows, 



dp 
di 



d_ 
dx 



d 5E d 5E\ 

1 dx 5p ^ m dx 5m J 



dm 
~dt 



d_ 

dx 



where we set 



m ii p 



Hp 



d_5^E 
} dx Sp 



and, as before, 



+ fj>f, 

dE 
dp 

m. 



d_SEY 
dx Sm J 



+ m x \i m x 



5E\ 



- 



5m J 



(9a) 
(9b) 



(1-«X 2 )" 



5E 
5m 



m. 



These equations have been introduced by Holm, Putkaradze and Tronci in [5], using a 
kinetic-theory description. The density and the magnetization vector are driven by the 
velocity 

d 5E d 5E 

V = H n — — + Hm. ■ — — . (10) 



^dx 5p +fim 



dx 5m 



The velocity advects the ratio \m\/p by 

d 



dt dx 



\m\ 
P 



We have the system energy 

E — \ I dxm ■ (l — a^dl) 1 m — | / dxp (l — a 2 p d^) 1 
Jn Jn 



(11) 



and, given a non-negative density, the second term is always non-positive. This represents 
an energy of attraction, and we therefore expect singularities in the magnetization vector to 
arise from a collapse of the particle density due to the ever-decreasing energy of attraction. 
There are three length scales in the problem that control the time evolution: the ranges a m 
and a p of the potentials in Eq. (11), and the smoothening length f3 m . 



Linear stability analysis 

We study the linear stability of the constant state (m, p) = (m Q ,p ). We evaluate the 
smoothened values of this constant solution as follows, 



d 2 f 

po = f(x)-a p ^, 
f (x) = po + Asinh (x/a p ) + B cosh (x/a p ) . 

For periodic or infinite boundary conditions, the constants A and B are in fact zero and 
thus 

(l-a p d 2 x y 1 p = p , (12) 
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and similarly po = (SE/Sm) mo = mo- The result (12) guarantees that the constant state 
(m ,po) is indeed a solution of Eq. Q. 

We study a solution (m, p) = (m + Sm, p + Sp), which represents a perturbation away 
from the constant state. By assuming that Sm and Sp are initially small in magnitude, we 
obtain the following linearized equations for the perturbation density and magnetization, 

d d 2 -i d 2 -i 



d_ 
dt 



5m = mo 



° 2 >- ,.2a2\-i e , d 2 , 2 



+ m x < m x 



9a ,2 { l - a l d D 5 P+ a^{ l - a mdl) rn -5m 

(1 - aldl)' 1 6m - (1 - Z^)" 1 5m] }. 



For mo 7^ we may choose two unit vectors hi and n 2 such that mo/\mo\, hi and 712 form an 
orthonormal triad (that is, we have effected a change of basis). We then study the quantities 
Sp, Sx, <5£i and 5£ 2 , where 



Sx = mo ■ 5m 
We obtain the linear equations 




££1 = ft-i • Sm, 5^2 = ^2 ■ Sw- 



Sp 



d 2 d 2 



0t l 



<9x 2 
5 2 



5X = -l^ol 2 ^ (1 -a'^ 2 ) X 5p+ \m \ 2 ^ (l-a 2 m d 2 x ) l S X , 



d 2 



.2 q2\"1 



<9x 



^=\{l-(3 2 m d 2 x y 1 -{l-a 2 m d 2 x y 1 } 5& 



1,2. 



By focusing on a single-mode disturbance with wave number we obtain the following 
system of equations 





(s P \ 
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dt 
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with eigenvalues 

CTq = 0, 



0"! 



p £; 2 



|m | 2 /c 2 



1 



1 



l + a 2 P 1 + c^A; 2 ' 



0"2 



1 + /^P 1 + a 2 m k 2 ■ 



(14) 



The eigenvalues are the growth rate of the disturbance (Sp, Sx, ££2) [22J. There are two 
routes to instability, when <7i > 0, or when o<i > 0. The first route leads to an instability 
when 

Po 



ui > 0, 



\mo\ 



> 



1 + a 2 p k 2 



l + al fc 2 ' 
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while the second route leads to instability when 

cr 2 > 0, a m > p m . 

We have plotted the growth rates for the case when p — l m o| 2 — 1> an d compared 
the theory with numerical simulations. There is excellent agreement at low wave numbers, 
although the numerical simulations become less accurate at high wave numbers. This can 
be remedied by increasing the resolution of the simulations. These plots are shown in 
Figs. [5] and [6j The growth rates o\p, are parabolic in k at small k; a± saturates at large 
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FIG. 5: (Color online) The first route to instability. Subfigure (a) shows the growth rate a\ for 
ot m < f3 m < OL p ^ with negativity indicating a stable equilibrium; (b) gives the growth rate G\ for 
ol p < a m < (3 m , with positivity indicating an unstable equilibrium. We have set |mo| = po = 1- 

k, while o"2 attains a maximum and decays at large k. The growth rates can be positive 
or negative, depending on the initial configuration, and on the relationship between the 
problem length scales. In contrast to some standard instabilities of pattern formation (e.g. 
Cahn-Hilliard [23] or Swift-Hohenberg [21]), the <Ti-unstable state becomes more unstable 
at higher wave numbers (smaller scales), thus preventing the 'freezing-out' of the instability 
by a reduction of the box size [23J. The growth at small scales is limited, however, by the 
saturation in a as k — > 00. Heuristically, this can be explained as follows: at higher wave 
number, the disturbance (Sp, Sx, 6£i, S&) gives rise to more and more peaks per unit length. 
This makes merging events increasingly likely, so that peaks combine to form larger peaks, 
enhancing the growth of the disturbance. 

Recall in Sec. [XX] that the different behaviors of the magnetization equation Q are the 
result of a competition between the length scales a m and (3 m . For a m < (3 m the initial (large- 
amplitude) disturbance tends to a constant, while for a m > (3 m the initial disturbance 
develops finer and finer scales. In this section, we have shown that the coupled density- 
magnetization equations are linearly stable when a m < /3 m , while the reverse case is unstable. 
In contrast to the first route to instability, the growth rate 02, if positive, admits a maximum. 
This is obtained by setting a' 2 (k) = 0. Then the maximum growth rate occurs at a scale 

Amax := 2nk m l x = 27Ta/ a m (5 m . 

Thus, the scale at which the disturbance is most unstable is determined by the geometric 
mean of a m and (3 m . Given a disturbance (<5p, 5x, 5£ 2 ) with a range of modes initially 
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FIG. 6: (Color online) The second route to instability. Subfigure (a) shows the growth rate o-i 
for a m < j3 m , with negativity indicating a stable equilibrium; (b) gives the growth rate a 2 for 
«m > Pm, with positivity indicating an unstable equilibrium. We have set \mo\ = po = 1. 



present, the instability selects the disturbance on the scale A max . This disturbance develops 
a large amplitude and a singular solution subsequently emerges. It is to this aspect of the 
problem that we now turn. 



Singular solutions 



In this section we show that a finite weighted sum of delta functions satisfies the par- 
tial differential equations (j9|. Each delta function has the interpretation of a particle or 
clumpon, whose weights and positions satisfy a finite set of ordinary differential equations. 
We investigate the two-clumpon case analytically and show that the clumpons tend to a 
state in which they merge, diverge, or are separated by a fixed distance. In each case, we 
determine the final state of the clumpon magnetization. 

To verify that singular solutions are possible, let us substitute the ansatz 



M 



M 



p (x, t) =} j a i (t) 8(x-Xj (t)) , m(x,t) = y ] bj (t) 8 (x - Xj (t)) , (15) 
i=i i=i 

into the weak form of equations ([9]). Here we sum over the different components of the 
singular solution (which we call clumpons). In this section we work on the infinite domain 
x E (—00,00). The weak form of the equations is obtained by testing Eqs. ^ with once- 
differentiable functions (x) and 1/; (x), 



d 
dt 



dxp (x, t) <p (x) 



, d SE 8 5E\ 

dx<P ( M ) (,,,__ + ^—^-1 



;i6a) 



d [ , / \ / / \ f 7 / / / \ , x ( d 5E d 6E\ 

— dxm (x, t) ■ ip (x) = - / dx^) [x) -m{x,t) \p p —— + p m - —— ) 



+ 



dxip (x) 



5E\ 

mx \ p m x — J 



(16b) 
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Substitution of the ansatz (15) into the weak equations (16) yields the relations 

ie{l,...,M}, (17) 



dcii 
~dt 



0. 



dt 



-V ( Xi 



dbj 
~dt 



biX [/IX 



6E\ 



8m J 



[X.-, 



where V and (fi x (5E/5m)) are obtained from the ansatz (15) and are evaluated at x^. 
Note that the density weights and the magnitude of the weights 6, remain constant in 
time. 

We develop further understanding of the clumpon dynamics by studying the two-clumpon 



version of Eqs. (17). Since the weights a%, a 2 , \b±\, and j^l are constant, two variables suffice 
to describe the interaction: the relative separation x = x± — x 2 of the clumpons, and the 
cosine of the angle between the clumpon magnetizations, cos<^ = b\ • &2 / [ &i 1 1 foa [ - Using the 
properties of the kernel E (0) = 1, E (0) = 0, we derive the equations 



dx 
~dt 



MH' ap (x) - B 1 H' am (x) H Pm (x) - B 2 H' am (x) y, y = cos <p 



dy 
dt 

where M = a\ + a 2 , B\ = |&i| 2 + \b 2 \ 2 , and B 2 = 2|&i||6 2 | are constants. Equations (18) form 
a dynamical system whose properties we now investigate using phase-plane analysis [25J. 
We note first of all that the \y\ > 1 region of the phase plane is forbidden, since the y- 



B 2 (l-y 2 ) [H Pm {x)-H am (x)] 
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FIG. 7: (Color online) The nullcline dx/dt = of the two-clumpon dynamical system with a m < a p . 
The region contained inside the dotted lines y = ±1 gives the allowed values of the dynamical 
variables (x,y). Subfigure (a) shows the case when j3 m < a m . The stable equilibria of the system 
are (x,y) = (±d, 1) and the line x = 0. All initial conditions flow into one of these equilibrium 
states; subfigure (b) shows the case when a m < (3 m . Initial conditions confined to the line y = 1 
flow into the fixed point (±c£, 1), while all other initial conditions flow into the line x = 0. 

component of the vector field (dx/dt, dy/dt) vanishes at \y\ = 1. The vertical lines x = and 
x = ±oo are equilibria, although their stability will depend on the value of the parameters 
(a m , a p , j3 m , B\, B 2 , M ). The curve across which dx/dt changes sign is called the nullcline. 



This is given by 



ME' 



B\E n 



X 



En 



B 2 E' 



(x) 



12 



4 










3 










2 














x= 


d/ 


\x=d 


1 







— f -* 


♦A- 



-1 










— >■ 




► 







< J.+. 


\x=d 

<-v£--->- — 








< — 








(a) 



(b) 



FIG. 8: (Color online) The nullcline dx/dt = of the two-clumpon dynamical system with a p < a m . 
The region contained inside the dotted lines y = ±1 gives the allowed values of the dynamical 
variables (x,y). Subfigure (a) shows the case when /3 m < a m . The lines x = and x = ±oo form 
the stable equilibria of the system. All initial conditions flow into one of these states; subfigure 
(b) shows the case when a m < (5 m . Initial conditions confined to the line y = 1 flow into the fixed 
points (0, 1) and (±oo, 1), while all other initial conditions flow into the line x = 0. 



which on the domain x G (— oo, oo) takes the form 
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-\x\ 



Bl -7T-M 
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Several qualitatively different behaviors are possible, depending on the magnitude of the val- 
ues taken by the parameters (a m , a p , /3 m , B x , B 2 , M). Here we outline four of these behavior 
types. 

• Case 1: The length scales are in the relation a m < Op, and (3 m < ot m . The vector field 
(dx/dt, dy/dt) and the nullcline are shown in Fig. [71(a). There is flow into the fixed 
points (x,y) = (±d, 1), and into the line x = 0, while y is a non-decreasing function of 
time, which follows from Eq. (18b). The ultimate state of the system is thus x = ±d, 
if = (alignment), or x = (merging). In the latter case the final orientation is given 



by the integral of Eq. (18b), 



tan 



tan 



2 



exp 



B 2 dt[H Pm ( X (t))-H am (x(t))\ 



fo 



<p(t = 0) 
(19) 



Case 2: The length scales are in the relation a m < a p , a m < f3 m . The vector field and 
the nullcline are shown in Fig. [7] (b). All flow not confined to the line y = 1 is into the 
line x = 0, since y is now a non-increasing function of time. The ultimate state of the 
system is thus x = ±d, tp — (alignment), or x = (merging). In the latter case the 
final orientation is given by the formula (19). 



Case 3: The length scales are in the relation a p < a m and (3 m < a m . The vector field 



and the nullcline are shown in Fig. |8| (a). Inside the region bounded by the line y = 
and the nullcline, the flow is into the line x = (merging), and the fixed points (±d, 1) 
are unstable. The flow below the line y — is towards the line x = 0. Outside of these 
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regions, however, the flow is into the lines x = ±00, which shows that for a suitable 
choice of parameters and initial conditions, the clumpons can be made to diverge. 

Case 4-' The length scales are in the relation ot p < a m and a m < (3 m . The vector field 



and the nullcline are shown in Fig. [8] (b). The quantity y is a non-increasing function 
of time. All flow along the line y = 1 is directed away from the fixed points (±d, 1) 
and is into the fixed points (0, 1), or (±00, 1). All other initial conditions flow into 
x = 0, although initial conditions that start above the curve formed by the nullcline 
flow in an arc and eventually reach a fixed point (x = 0, y < 0). 

We summarize the cases we have discussed in Table |TT| Using numerical simulations of 
Eqs. (18), we have verified that Cases (1)— (4) do indeed occur. The list of cases we have 
considered is not exhaustive: depending on the parameters B\, B 2 , and M, other phase 
portraits may arise. Indeed, it is clear from Fig. [7] that through saddle-node bifurcations, 
the fixed points (x,y) = (±d, 1) may disappear, or additional fixed points (x,y) = (±<f, —1) 
may appear. Our analysis shows, however, that it is possible to choose a set of parameters 



Case 


a m vs. a p 


OL m vs. m 


Equilibria 


Flow 


(1) 
(2) 
(3) 
(4) 


Om ^ Otp 

dp < a m 


0m ^ O m 
*^ 0m 

0m <~ 0! m 
^ 0m 


(x, y) = (±d, 1); x = 0; x = ±00 
(x, y) = (±d, 1); x = 0; x = ±00 
(x, y) = (±d, 1); x = 0; x = ±00 
(x, y) = (±d, 1); x = 0; x = ±00 


Flow into x = and (x, y) = (±d, 1) 
Flow into x = and (x, y) = (±d, 1) 
Flow into x = and x = ±00 
Flow into x = and x = ±00 



TABLE II: Summary of the distinct phase portraits of Eq. (18) studied. 



i ra , m , Bi, B 2 , M) such that two clumpons either merge, diverge, or are separated by 
a fixed distance. 



Numerical Simulations 



To examine the emergence and subsequent interaction of the clumpons, we carry out 
numerical simulations of Eq. ^ for a variety of initial conditions. We use an explicit finite- 
difference algorithm with a small amount of artifical diffusion. We solve the following weak 
form of Eq. ([9]), obtained by testing Eq. (16) with Hp m , 



d 2 p f 

Atrtif7^+ / dyH Pm (x-y)p(y,t)V(y,t) 



dp n d 2 p 
di 



d^j _ n d 2 pi 
dt artif dx 2 + 



+ / d y H Pm i x - y) e i 



SE\ 



where p = Hp m * p and is the unit vector in the i th direction. We work on a periodic 
domain Q = [— L/2, L/2], at a resolution of 250 gridpoints; going to higher resolution does 
not noticeably increase the accuracy of the results. 
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The first set of initial conditions we study is the following, 

m (x, 0) = (sin (4k x + (p x ) , sin (4/c x + <f> y ) , sin (4k x + 2 )) , 
p(x,0) = 0.5 + 0.35 cos (2kox), 



(20) 



where <f> x , <f> y , and <p z are random phases in the interval [0, 27r], and k = 2ir/L is the 
fundamental wave number. The initial conditions for the magnetization vector are chosen to 
represent the lack of a preferred direction in the problem. The time evolution of equations ffify 
for this set of initial conditions is shown in Fig. [9j After a short time, the initial data become 
singular, and subsequently, the solution (p, m) can be represented as a sum of clumpons, 



M 



M 



p (x, t) = di5 (x - Xi (t)) , m (x, t) = 6j (t) 5 (rr - (£)) 



M 



i=l 



Here M = 2 is the number of clumpons present at the singularity time. This number 
corresponds to the number of maxima in the initial density profile. The forces exerted by 
each clumpon on the other balance because of the effect of the periodic boundary conditions. 
Indeed, any number of equally-spaced, identical, interacting particles arranged on a ring are 
in equilibrium, although this equilibrium is unstable for an attractive force. Thus, at late 
times, the clumpons are stationary, while the magnetization vector p shows alignment of 
clumpon magnetizations. 
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FIG. 9: (Color online) Evolution of sinusoidally-varying initial conditions for the density and 



magnetization, as in Eq. (20). Subfigure (a) shows the evolution of H p * p for t £ [0, 0.15], by which 



time the initial data have formed two clumpons; (b) shows the evolution of p x . The profiles of p y 
and p z are similar. Note that the peaks in the density profile correspond to the troughs in the 
magnetization profile. This agrees with the linear stability analysis, wherein disturbances in the 
density give rise to disturbances in the magnetization. 

We gain further understanding of the formation of singular solutions by studying the 
system velocity V just before the onset of the singularity. This is done in Fig. 10 Fig- 
ure 



10 



(a) shows the development of the two clumpons from the initial data. Across each 
density maximum, the velocity has the profile V ~ A (t) x, where A (t) > is an increasing 
function of time. This calls to mind the advection problem for the scalar 9 (x,t), studied by 
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(a) 



(b) 



FIG. 10: (Color online) Evolution of sinusoidally-varying initial conditions for the density and 



magnetization, as in Eq. (20). Subfigure (a) shows the system velocity V given in Eq. (10), just 
before the singularity time; (b) shows the magnetization p at the same time. The density maxima 
emerge at the locations where the convergence of — V (flow into x = and x = ±L/2) occurs, and 
the magnetization develops extrema there. 



Batchelor in the context of passive-scalar mixing 

de.de . n 

— = A x— , A > 0. 
ot ox 

Given initial data 9 (x, 0) = 9oe~ x2 ^ e o, the solution evolves in time as 

9(x,t) = 9 e~* 2 ^ e - 2X " t ), 
so that gradients are amplified exponentially in time, 

dx~ t\ 

in a similar manner to the problem studied. 



The evolution of the set of initial conditions (20 ) has therefore demonstrated the following: 
the local velocity V is such that before the onset of the singularity, matter is compressed 
into regions where p (x, 0) is large, to such an extent that the matter eventually accumulates 
at isolated points, and the singular solution emerges. Moreover, the density maxima, rather 
than the magnetization extrema, drive the formation of singularities. This is not surprising, 
given that the attractive part of the system's energy comes from density variations. 

To highlight the interaction between clumpons, we examine the following set of initial 
conditions, 

m (x, 0) = m = const., p (x, 0) = 0.5 + 0.35 cos (Sk x) , (21) 

where ko = 2n / L is the fundamental wave number. Since this set of initial conditions 
contains a large number of density maxima, we expect a large number of closely-spaced 
clumpons to emerge, and this will illuminate the clumpon interactions. The time evolution 
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FIG. 11: (Color online) Evolution of a fiat magnetization field and a sinusoidally- varying density, 



as in Eq. ( |21[ ). Subfigure (a) shows the evolution of H p * p for t G [0.5, 1]; (b) shows the evolution 
of \x x . The profiles of \i y and fj, z are similar. At t = 0.5, the initial data have formed eight equally 
spaced, identical clumpons, corresponding to the eight density maxima in the initial configuration. 
By impulsively shifting the clumpon at x = by a small amount, the equilibrium is disrupted and 
the clumpons merge repeatedly until only one clumpon remains. 



of equations ^ for this set of initial conditions is shown in Fig. 11 As before, the solution 
becomes singular after a short time, and is subsequently represented by a sum of clumpons, 



M 



p (x, t) = aid (x — Xi (t)) , m (x, t) 



i=i 



M 

£ 

4=1 



h (t)6(x- Xi (t)) 



M = 8. 



Here M = 8 is the number of clumpons at the singularity time. This number corresponds to 
the number of maxima in the initial density profile. As before, this configuration of equally 
spaced, identical clumpons is an equilibrium state, due to periodic boundary conditions. 
Therefore, once the particle-like state has formed, we impulsively shift the clumpon at 
x = by a small amount, and precipitate the merging of clumpons. The eight clumpons 
then merge repeatedly until only a single clumpon remains. The tendency for the clumpons 
to merge is explained by the velocity V, which changes sign across a clumpon. Thus, if 
a clumpon is within the range of the force exerted by its neighbours, the local velocity, if 
unbalanced, will advect a given clumpon in the direction of one of its neighbours, and the 



clumpons merge. This process is shown in Fig. 13 



IV. CONCLUSIONS 

We have investigated the non-local Gilbert (NG) equation introduced by Holm, 
Putkaradze, and Tronci in [6| using a combination of numerical simulations and simple an- 
alytical arguments. The NG equation contains two competing length scales of non-locality: 
there is a length scale a associated with the range of the interaction potential, and a length 
scale (3 that governs the smoothened magnetization vector that appears in the equation. 
When a < (5 all initial configurations of the magnetization tend to a constant value; while 
for (3 < a the initial configuration of the magnetization field develops finer and finer scales. 
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FIG. 12: (Color online) Evolution of a flat magnetization field and a sinusoidally-varying density, as 



in Eq. (21 ). Subfigure (a) shows the system velocity V given in Eq. (10) just before the singularity 
time; (b) gives the magnetization [i at the same time. The density maxima emerge at the locations 
where the convergence of — V occurs, and the magnetization develops extrema there. 
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FIG. 13: (Color online) Evolution of a flat magnetization field and a sinusoidally-varying density, 



as in Eq. (21 ). Shown is the velocity profile for t £ [0.5, 1]; the system velocity is given by Eq. (10). 
The velocity —V flows into each density maximum, concentrating matter at isolated points and 
precipitating the formation of eight equally-spaced identical clumpons. On a periodic domain, such 
an arrangement is an equilibrium state, although it is unstable. Thus, by impulsively shifting the 
clumpon at x = by a small amount, we force the clumpons to collapse into larger clumpons, until 
only a single clumpon remains. 



These two effects are in balance when a = (3, and the system does not evolve away from 
its initial state. Furthermore, the NG equation conserves the norm of the magnetization 
vector m, thus providing a pointwise bound on the solution and preventing the formation 
of singular solutions. 

To study the formation of singular solutions, we couple the NG equation to a scalar 
density equation. Associated with the scalar density is a negative energy of attraction that 
drives the formation of singular solutions and breaks the pointwise bound on the m. Three 
length scales of non-locality now enter into the problem: the range of the force associated 
with the scalar density, the range of the force due to the magnetization, and the smoothening 
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length. As before, the competition of length scales is crucial to the evolution of the system; 
this is seen in the linear stability analysis of the coupled equations, in which the relative 
magnitude of the length scales determines the stability or otherwise of a constant state. 

Using numerical simulations, we have demonstrated the emergence of singular solutions 
from smooth initial data, and have explained this behavior by the negative energy of attrac- 
tion produced by the scalar density. The singular solution consists of a weighted sum of delta 
functions, given in Eq. (15), which we interpret as interacting particles or clumpons. The 
clumpons evolve under simple finite-dimensional dynamics. We have shown that a system 
of two clumpons is governed by a two-dimensional dynamical system that has a multiplicity 
of steady states. Depending on the length scales of non-locality and the clumpon weights, 
the two clumpons can merge, diverge, or align and remain separated by a fixed distance. 

Our paper thus gives a qualitative description of the dynamics. Future work will focus on 
the regularity of solutions of the NG equation, and the existence and regularity of solutions 
for the coupled density-magnetization equations. Bertozzi and Laurent [27] have studied 
the simpler (uncoupled) non-local scalar density equation, proving existence, uniqueness, 
and blowup results using techniques from functional analysis, and a similar analysis will 
illuminate the equations we have studied. The behavior of singular solutions in higher 
dimensions is another topic that deserves further study. 
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